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ABSTRACT 


An  exciting  prospect  for  modern  filtering  theory  lies  in  the  field 
of  orbital  navigation,  bat  as  of  today  the  most  impressive  results 
in  orbit  determination  have  been  obtained  by  classical  methods  of 
parameter  estimation.  Some  current  experience  and  problems 
are  described  and  their  implications  for  modern  filtering  theory 
are  discussed.  Two  techniques,  found  effective  in  the  classical 
process,  are  presented  in  the  context  of  recursive  estimation. 
Four  illustrative  plots  are  included. 
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1.  INTRODUCTION 


Many  early  expositions  of  the  modern  filtering  theory  (MFT)  have  been  in 
orbit  determination  and  some  of  the  most  exciting  prospects  for  the  imple¬ 
mentation  of  MFT  lie  in  the  field  of  space  navigation  (Refs.  1  and  2).  How¬ 
ever,  as  of  today,  the  most  impressive  results  in  orbit  determination  have 
been  obtained  with  the  classical  methods  of  parameter  estimation,  as  applied 
at  the  Naval  Weapons  Laboratory,  Applied  Physics  Laboratory  (APL)  of  The 
Johns  Hopkins  University,  and  at  other  installations  (Refs.  3,  4,  and  5). 

Some  of  the  experience  and  problems  of  the  classical  approach  will  be  de¬ 
scribed  and  their  implications  for  MFT  discussed.  Two  techniques  used 
effectively  in  the  classical  process  are  presented  in  the  context  of  recursive 
estimation. 
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2.  THE  PROBLEMS 


The  typical  problems  of  MFT  (smoothing,  filtering,  and  prediction)  as  defined 
by  Kalman  (Ref.  6)  correspond  in  the  field  of  orbit  determination  to  those  of 
ephemeris  reconstruction,  navigation,  and  orbital  prediction,  all  of  -which 
have  important  applications.  The  navigational  or  filtering  problem  (that  of 
seeking  a  current  estimate  of  position  and  velocity  based  on  a  current  obser¬ 
vation  and  a  recent  past  estimate)  is  one  that  astronauts  will  face  with  the 
advent  of  onboard  navigation,  but  it  has  not  commanded  much  specific  atten¬ 
tion  in  classical  orbit  determination.  As  for  the  smoothing  and  prediction 
problems,  the  ephemeris  of  a  geodetic  satellite  must  be  reconstructed  with, 
all  possible  precision,  and  a  navigational  satellite  is  useful  only  if  its  position 
can  be  predicted  with  great  accuracy.  There  have  been  impressive  results  • 
in  both  areas. 
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3.  THE  RECORD 


In  both  problems  described,  primary  concern  is  with  the  satellite  ephemeris, 
that  is,  with  its  position  and  velocity  as  Sanctions  of  time.  A  variety  of 
empirical  methods  have  been  used  to  increase  the  accuracy  of  the  ephemeris. 
Nearly  all  were  found  wanting  and  replaced,  with  improved  results,  by  more 
sophisticated  estimation  procedures. 


Following  are  two  related  examples:  In  some  early  attempts  at  orbit  recon¬ 
struction  at  Aerospace  Corp.  ,  the  final  residuals  from  several  flights  were 
much  too  large  (several  thousands  of  feet)  for  the  quality  of  the  data,  but 
could  be  accounted  for  very  nicely  with  a  time  correction.  Timing  error 
curves  were  prepared  and  used  as  an  empirical  device  for  increasing  the 
accuracy  of  the  ephemeris.  Their  shape  (sinusoidal  with  a  12-hour  period) 
strongly  suggested  that  the  principal  longitude-dependent  term  (the  J22  ‘firm) 
should  be  included  in  the  geopotential.  This  was  done ,  and  a  suitable  value  for 
J22  did  indeed  effect  a  dramatic  reduction  in  the  residuals.  The  acid  test  of 
a  scientific  hypothesis  being  prediction,  the  revised  geopotential,  if  a  true 
scientific  advance  and  not  just  a  fancier  empirical  device,  would  have  to  lead 
to  smaller  residuals  on  other  flights  as  well.  It  did,  somewhat,  but  for  best 
results  required  a  different  value  for  J22,  again  due  to  further  deficiencies 
in  the  geopotential  model.  We  were  soon  in  the  business  of  estimating  J22 
and  a  few  other  coefficients  in  the  geopotential,  as  well  as  a  drag  coefficient, 
in  addition  to  the  six  orbit  parameters.  All  this  was  in  quest  of  just  a  position 
and  velocity  ephemeris. 


Additional  degrees  of  freedom  will  invariably  reduce  the  residuals  in  the  fitted 
interval,  but  this  is  no  proof  that  a  better  model  has  been  found.  As  men¬ 
tioned,  the  proper  test  of  a  hypothesis  is  in  prediction;  and  the  use  of  empirical 
techniques  to  obtain  a  good  reconstruction  has  invariably  resulted  in  rapidly 
increasing  residuals  outside  the  fitted  interval.  Ultimately,  only  the  true 
model  will  pass  the  test  of  prediction. 
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The  outstanding  example  of  precision  orbit  determination  today  is,  possibly, 
the  work  on  the  TRANSIT  project  by  APLi.  The  project  requires  the  predic¬ 
tion  of  a  precision  position  and  velocity  ephemeris,  and  thus  a  highly  accurate 
knowledge  of  the  geopotential,  {Their  satellites  are  sufficiently  high  so  that 
drag,  although  it  cannot  be  ignored,  is  not  a  problem. )  To  fulfill  their  mis¬ 
sion,  they  have  had  the  huge  problem  of  estimating  70  non-zonal  parameters 
of  the  geopotential  with  1662  passes  of  satellite  data  (from  5  satellites)  taken 
at  uncertain  locations  (Refs.  4  and  5). 

Working  with  the  same  data,  but  with  independent  computer  programs,  the 
Naval  Weapons  Laboratory  has  also  made  an  impressive  determination  of  the 
coefficients  of  the  geopotential  (Ref.  3).  In  their  program,  as  many  as  500 
parameters  are  estimated  simultaneously. 

Inspection  of  residuals  played  an  important  role  in  their  work,  as  it  has  in 
ours.  Noting  a  pronounced  oscillation  of  period  about  60  hours,  APL  deduced 
a  resonance  effect  between  the  orbital  period  of  the  satellite  and  non-zonal 
harmonics  in  the  geopotential  of  order  13,  with  a  beat  period  of  60  hours 
(Ref.  7).  They  solved  for  the  appropriate  geopotential  coefficients,  and  noted 
an  appropriate  reduction  in  residuals  on  two  satellites  launched  three  months 
apart  (the  prediction  test),  thereby  completing  an  iteration  of  the  hypothesis- 
solution-test  cycle. 

The  remaining  residuals  exhibit  a  daily  pattern  correlated  with  the  longitude 
of  the  ascending  node  and,  for  some  satellites,  with  the  latitude.  Only  10  to 
30  m  of  the  remaining  residuals  (less  than  80  m)  appear  to  be  essentially 
random. 

These  achievements  are  very  impressive.  Not  only  are  good  fits  obtained  in 
reconstruction,  but  also  their  results  satisfy  several  important  and  indepen¬ 
dent  checks,  and  pass  the  test  of  prediction  brilliantly. 
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4.  ESTIMATING  UNCERTAINTIES 


There  are  two  types  of  measures  of  the  uncertainty  in  a  reconstructed 
ephemeris.  One  (which  might  be  called  "practical")  is  derived  from  an 
examination  of  the  residuals;  the  other,  the  covariance  matrix,  is  (under 

sjc 

certain  conditions)  the  inverse  of  the  normal  matrix  A'WA,  and  is  thus  a 
theoretical  measure  in  which  the  residuals  have  no  part. 

It  is  easy  to  stage  examples  wherein  small  changes  in  the  geophysical  model 
can  have  large  effects  upon  the  residuals,  but  virtually  none  upon  the  covari¬ 
ance  matrix.  These  do  not  speak  well  for  the  value  of  the  covariance  matrix, 
as  a  measure  of  uncertainty.  The  explanation  is  that  the  conditions  under 
which  the  covariance  matrix  is  valid  are  overlooked  -  as  usual.  The  conditions 
are  that  the  random  observational  errors  have  mean  zero  and  covariance 
matrix  W  ^  and  that  the  model  be  correct. 

The  latter  requirement  is  invariably  violated,  in  orbit  determination  problems, 
by  linearization  of  the  original  equations,  and  more  subtly,  every  time  a 
geopotential  or  atmospheric  model  is  used;  for,  inevitably,  neither  the  true 
forma  of  the  models  nor  the  coefficients  involved  are  certainly  known.  Nor 
are  the  locations  of  the  observing  stations  or  the  systematic  error  charac¬ 
teristics  of  the  observing  instruments  precisely  known.  Thus  the  usual 
covariance  matrix  determined  in  the  course  of  estimating  six  (or  a  few  more) 
orbit  parameters,  assuming  all  other  parameters  are  exactly  known,  may  be 
(depending  on  the  sophistication  and  accuracy  of  the  model)  hopelessly  opti¬ 
mistic  and  virtually  worthless  as  a  measure  of  ephemeris  accuracy. 

Yet  it  is  precisely  this  matrix  that  governs  the  size  of  corrections  to  MFT 
state  variables.  A  small  covariance  matrix,  implying  a  high  but  false  confi¬ 
dence  in  the  current  estimates,  will  permit  only  small  corrections  thereto. 
Clearly,  some  techniques  that  tend  to  de-weight  old  observations  will  be 
required. 


The  prime  denotes  the  transpose  of  a  vector  or  matrix. 
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For  system  design  studies,  however,  there  is  no  recourse  to  uncertainty 
estimates  derived  from  observations.  Two  alternatives  are  a  Monte  Carlo 
analysis  (of  frightening  dimensions  for  a  typical  orbit  determination  problem), 
or  an  expanded  form  of  covariance  analysis  described  in  Refs.  8,  9,  and  10 
and  below. 

It  is  assumed  here  that  the  form  of  a  linear  model  is  known,  and  that  certain 
parameters  or  state  variables  are  to  be  estimated,  but  in  the  presence  of 
other  parameters  the  values  of  which  are  uncertain  and  not  being  estimated. 
From  the  first  assumption,  it  is  apparent  that  the  method  has  some  deficien¬ 
cies,  but  may  permit  sensible  answers  to  relative  questions  such  as;-  "How 
does  the  uncertainty  in  the  range  bias  of  Station  A,  which  is  not  being  esti¬ 
mated,  degrade  our  positional  accuracy?" 

Suppose  that  observational  residuals  can  be  represented,  after  linearization, 
as  the  vector 


y  =  Ax  +  Bv  +  c 


(1) 


where 

y  is  an  n-vector  of  residuals 

x  is  a  p-vector  of  parameters  to  be  estimated 

v  is  the  q-vector  of  parameters  not  being  estimated  but 
with  uncertainties  given  by  E(w' )  - 

A,B  are  matrices  (n  X  p  and  n  X  q,  respectively) 

€  is  an  n-vector  of  random  observational  errors,  of  zero 
mean  and  covariance  matrix  E(t  c ')  =  W-* 

J 

In  deriving  the  estimate  x  of  x,  however,  the  dependence  upon  v  is  ignored. 
Thus  the  estimate  is  the  familiar 

x  =  (A'WAJ^A'Wy  (2) 
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which  is  readily  seen  to  be  an  unbiased  estimate  of  x  if  E(v)  =  0  (meaning 
that  unbiased  estimates  were  used  as  nominal  values  of  the  unknown  param¬ 
eters).  Then 

5x  =  x  -  x  =(A'WA)‘1A'W(Bv+  t)  (3) 

and  the  covariance  matrix  of  the  estimate  is,  assuming  5(vt ')  =  0, 

E(6xox')  =  (A'WA)*1  +  (A'  WA)_1(A'  WB)E(w'  ){A'  WB)'  (A'  WA)"1  ,  (4) 


or 


P 

XV 


P  +PQP  Q'P 

X  X  V  X 


Now  P^.  and  Q  are  readily  formed  recursively.  The  matrix  Px  =  (A'WA) 
can  be  computed  recursively,  using  an  identity,  without  explicit  matrix 
inversion  if  W  is  diagonal  (Refs.  2  or  8): 


(5) 

-1 


p*<k  +  »  -  p,<ki  -  *  ’iii]'1  w*<*>  <« 

Also 


Q(k  +  1}  -  Q(k)  +  ak+lWk+lbk+l 


(7) 


where 


Px(k),  Q(k) 


“k+l*  ^+1 
-1 

wk+l 


signifies  the  evaluation  of  the  matrices  using  the  first  k 
rows  of  A,  B,  W 

are  the  (k  +  l)8t  rows  of  A,  3 

is  the  (k  +  l8*  diagonal  element  of  W  * 
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We  fee  that  the  uncertainty  inx,  as  given  by  P  ,  is  a  sum  of  two  terms,  the 
second  arising  from  the  use  of  uncertain  values  for  v.  The  second  term  can 
be  further  expanded  by  partitioning  the  matrices  involved,  to  present  sepa¬ 
rately  the  effects  upon  of  uncertainties  in  individual,  or  groups  of, 
parameters.  Assume  that  P^  is  block  diagonal 


P  =|P  .,  P  ,, 

v  I  v ,  1  *  v,  2 


r,  sj 


(8) 


and  specify  a  compatible  partition  of  B, 


B  =  [Blt  B2,  •  •  -  ,  Bg] 
and  let  Q.  =  A'WB.(i  =  1,  •  *  •  ,  s).  Then 


(9) 


P  =P  +P[Q.P  .q;  +  •  •  •  +  Q  P  Q']P 

XV  X  X  1  V,  1  1  S  V,  S  S  X 


(10) 


and  P  Q.P  .Q.'P  is  the  contribution  to  the  covariance  matrix  P  due  to 

X  1  V,  1  1  X  ^  3CV 

the  uncertainties  in  the  i  group  of  non-estimated  parameters. 
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5.  DATA  EDITING  AND  BOUNDED  ESTIMATES 


Another  real  problem,  which  will  have  to  be  faced  in  orbital  applications  of 
MFT,  is  that  of  data  editing.  In  the  techniques  used  with  the  classical  methods 
of  orbit  determination  it  is  usually  assumed  that  many  observations  are 
simultaneously  available  and  that,  by  polynomial  or  other  smoothing,  outliers 
can  be  detected  and  discarded.  The  assumption  is,  however,  contrary  to  the 
spirit  of  MFT.  If  observations  are  comparatively  scarce  and  only  available 
singly,  the  usual  kind  of  data  editing  may  be  out  of  the  question.  There  is  the 
further  possibility  that  an  apparently  bad  observation  is  in  fact  valid,  the 
large  residual  being  due  to  the  inadequacy  of  the  mathematical  model.  In  such 
circumstances,  a  practical  procedure  might  be  to  reject  those  observations 
that  induce  an  unacceptably  large  correction  to  the  estimate.  An  alternative 
would  be  to  accept  all  observations,  but  to  impose  a  bound  upon  the  computed 
correction  to  the  estimate. 

th 

Suppose  that,  after  recording  the  k  observation,  we  seek  a  correction 
Ax(k)  =  x(k)  -  x(k  -  1)  to  the  parameter  vector  that  minimizes 

f(k)  =  [y(k)  -  A(k)x(k)]'W(k)[y(k)  -  A(k)x(k)]  (11} 

subject  to  the  bound 


|G*Ax(k)  I 


■id)- 


Ax(k) '  G^Ax(k)  5  1 


(12) 


i=l 


-1 


(G  is  a  diagonal  matrix  with  elements  . ) 

The  solution  of  the  bounded  problem  requires  the  minimization  of 

f*(k)  =  f(k)  +  \2(k)Ax(k)'G2Ax(k) 


(13) 
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where  \2{k)  is  a  positive  constant  to  be  determined.  (The  method  is  due  to 
D.  Morrison  and  reported  in  Refs.  10  and  11.)  To  locate  the  minimum  of 
f*(k),  we  set  3f #{k) / 3 [  Ax(k)]  =  0  and  solve  (with  \(k)2  =  0  initially}  the  result¬ 
ing  linear  system 

[A(k)'W(k)A(k)  +  X2(k)G2]Ax(k) 

=  A(k)'W(k)[y(k)  -  A(k)x(k  -  1)]  (14) 

■  t**  -  ‘>'w<k  -  --1’] 

(the  dotted  bars  indicate  a  partitioning) 

=  A(k  -  1)' W(k  -  l)[y(k  -  1)  -  AU  -  l)x(k  -  I)]  +  a£wk[yk  -  a^fr  -  1)] 

(15) 

=  a^wk[yk  -  a^Ck  -  1)]  (16) 

Since  the  first  member  of  Eq.  (15)  will  be  zero  as  the  result  of  minimizing 

f(k  -  1). 

If  1h.e  solution  fails  to  satisfy  the  bounding  condition,  a  non-zero  value  of 
2 

X  (k)  is  chosen  to  initiate  a  search  and  iteration  procedure  for  a  value  which 

z 

does  satisfy  |[ G  •  Ax(k)  ||  =  1,  at  least  approximately.  Each  step  requires  tire 

solution  of  the  linear  system,  and,  upon  convergence,  the  solution  minimizes 
* 

f  (k)  subject  to  the  bounding  condition  in  Sq.  (12). 

5. 1  APPROXIMATE  INVERSION 

Recall  that  P  =  (A'WA)"1  can  be  found  recursively  through  the  use  of  a  matrix 
identity.  By  another  application  Of  the  same  formula,  we  obtain 

[P-1  +  (XG)'I(XG)]-1  =  P  -  X2FG(X2GPG  +  I)-1GP  (17) 
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in  which  the  inverse  can  be  approximated  by  truncating  a  series  representation 


(I  +  X2GPG)_1  =  I  -  \2GPG  +  (\2GPG)2  -  .  .  • 


(18) 


leading  to 


[P'1  +  X2G2]-1  =  P[I  -  C  +  C2  -  .  .  .  ]  (19) 

2  2 

where  C  =  X  G  P.  Thus,  even  the  bounded  problem  may  permit  a  solution 
(although  only  approximate)  without  explicit  matrix  inversion. 

5.2  SUCCESSIVE  BOUNDING 

Bounds  may  he  imposed  in  successive  stages,  in  which  case  the  right  member 
of  Eq.  (lb)  requires  modification.  Had  the  estimate  x(k  -  1)  been  the  result 

3*£ 

of  minimizing  f  (k  -  1),  [that  is,  had  Ax(k  -  1)  also  been  the  solution  of  an 
equation  like  Eq.  (14)],  then  in  Eq.  (15)  for  the  stage 

A(k  -  l)'W(k  -  l)[y(k  -  1)  -  A(k  -  l)x(k  -  1)]  =  \2(k  -  l)G2Ax{k  -  1)  (20) 

and  this  quantity  must  be  added  to  Eq.  (16),  or  finally 

[A(k)'W(k)A(k)  +  X2(k)G2]Ax(k)  =  a^wk[yk  -  ajX(k  -  1)]  +  X2(k  -  l)G2dx(k  -  1) 

(21) 


-13- 


6.  SUMMARY  .AND  SPECULATION 


Orbit  determination  accuracies  have  improved  dramatically  in  eight  years  - 
from  the  early  days  of  estimating  only  orbit  parameters,  with  simple  models 
and  empirical  techniques  -  to  today's  sophisticated  methods  wherein  realistic 
complex  models  can  be  used  and  their  many  parameters  estimated.  The 
emphasis  has  been  on  analysis  of  residuals,  explanatory  deterministic 
hypotheses,  and  testing  by  prediction.  The  random  forcing  function,  common 
in  the  dynamic  problems  of  MFT ,  has  played  no  significant  part  in  modern 
orbit  determination  of  the  classical  type. 

The  drag  problem,  especially  that  of  determining  short -period  fluctuations 
in  density,  may  be  a  candidate  for  stochastic  treatment.  Rauch  (Ref.  12)  has 
obtained  nice  reductions  in  residuals  in  a  smoothing  application  of  MFT,  but 
it  is  hard  to  foresee  effective  prediction  from  this  approach.  Perhaps  the 
best  use  of  MFT  methods  will  he  in  the  estimation  of  amplitudes,  periods, 
and  time  constants  as  data  for  the  formulation  of  deterministic  hypotheses. 
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APPENDIX.  EXAMPLES 


All  the  results  displayed  were  obtained  using  the  program  TRACE  (Ref.  10), 
which  can.  generate  data,  perform  covariance  analyses,  and  estimate  trajec¬ 
tory,  geophysical,  and  observational  parameters.  The  estimation  methods 
are  basically  classical  and  operation  of  the  program  in  a  filtering  mode  is 
inefficient  and  not  feasible  for  extensive  studies. 

The  reference  trajectory  for  all  the  figures  is  a  circular  equatorial  orbit  at 
am  altitude  of  100  n  mi  above  a  spherical  earth.  The  only  perturbative  force 
is  that  of  atmospheric  drag  derived  from  a  Jacchia  model  atmosphere  and  a 
ballistic  coefficient  C^A/W  =  0.  1  ft  /lb.  Observing  stations  are  located  at 
30*  intervals  in  the  orbital  (equatorial)  plane.  By  specifying  an  interval  of 
two  minutes,  with  a  minimum  elevation  of  20*,  only  one  observation  per  pass 
was  generated,  to  which  independent  normal  random  noise  was  added. 

Two  revolutions  of  data  were  fitted  with  both  the  true  and  a  false  (C^A/W  =  0.  11) 
model.  The  in-track  separation  between  the  reference  and  the  fitted  trajec¬ 
tories  for  the  fitting  interval,  and  for  a  like  two -revolution  prediction  interval, 
is  shown  in  Fig.  1.  By  using  the  true  model,  the  maximum  amplitudes  in  the 
reconstruction  and  prediction  intervals  are  142  and  364  ft,  respectively;  the 
curve  would  lie  entirely  within  the  shaded  area. 

In  Fig.  2  the  residuals  for  a  false -model  fit  and  prediction  are  displayed. 

Note  their  systematic  character,  the  level  in  the  reconstruction  interval 
(10  times  the  random  component),  and  the  rapid  growth  in  the  prediction 
interval. 

When  the  model  parameters  are  assumed  known,  the  covariance  matrices 
give  an  in -track  standard  deviation  of  only  200  ft  at  the  end  of  the  prediction 
interval.  The  10-percent  change  in  C^A/W,  which  leads  to  a  9-mile  difference 
in  trajectory  position,  causes  only  a  4 -ft  change  in  the  in-track  standard 
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Figure  1.  In-Track  Differences 


nmi 


-£1- 


Figure  2.  Range  Residuals 


deviation.  Much  more  realistic  statistics  are  produced  when  a  10-percent 
uncertainty  is  introduced  for  C^A/W;  that  is,  Pv  =  0.  0001  in  Eq.  (5).  The 
results  are  shown  in  Fig.  3,  in  which  200  ft  are  represented  by  about  half 
the  thickness  of  the  shaded  area. 

Operating  TRACE  in  recursive  modes  led  to  the  trajectory  differences  (from 
the  reference  trajectory)  displayed  in  Fig.  4.  The  theoretical  values  are 
taken  from  covariance  calculations  and,  in  the  prediction  interval,  are  simply 
those  of  Fig.  3  presented  in  logarithmic  scale.  Again,  the  trajectory  deter¬ 
mined  in  the  presence  of  a  false  (by  10  percent  in  C^A/W)  model  diverged 
rapidly  from  the  reference  trajectory,  but  the  differences  were  accurately 
predicted  by  the  covariance  matrix. 

Two  runs  wer^  made  in  which  the  correct  value  of  C^A/W  was  used,  and  the 
different  results  illustrate  the  effect  of  the  assumption  that  the  orbit  filtering 
problem  can  be  linearized.  In  the  linear  mode,  TRACE  is  mathematically  a 
static  model  recursive  estimation  program.  The  initial  position  and  velocity 
parameters  are  re-estimated  with  each  observation,  as  per  Eq.  (16),  with 
X^(k)  =0.  In  both  the  linear  and  nonlinear  modes,  a)  the  covariance  matrix 
Px  is  obtained  by  direction  inversion  rather  than  by  appeal  to  the  identity  in 
Eq.  (6),  and  b)  the  position  of  the  satellite  is  computed  from  the  full  (nonlinear) 
equations  of  motion.  In  the  nonlinear  mode,  TRACE  successively  fits  batches 
of  1,  2,  3  .  ,  .  observations  without  iteration.  The  matrix  A  and  the  residuals 
y  are  recomputed  for  each  fit  and  thus  the  nonlinearities  in  the  covariance 
matrix  and  residual  vector  are  taken  into  account.  Standard  deviations  of 
20,000  ft  and  40  fps  were  used  in  constructing  the  initial  covariance  matrix 
required  for  recursive  estimation. 
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igure  3.  In-Track  Standard  Deviation 
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Figure  4.  Filtered  Trajectory  Differences 
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